SIMULTANEOUS RECONSTRUCTION OF OUTER BOUNDARY 
SHAPE AND ADMITTIVITY DISTRIBUTION IN ELECTRICAL 
IMPEDANCE TOMOGRAPHY 
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Abstract. The aim of electrical impedance tomography is to reconstruct the admittivity dis- 
tribution inside a physical body from boundary measurements of current and voltage. Due to the 
severe ill-posedness of the underlying inverse problem, the functionality of impedance tomography 
relies heavily on accurate modelling of the measurement geometry. In particular, almost all recon- 
struction algorithms require the precise shape of the imaged body as an input. In this work, the 
need for prior geometric information is relaxed by introducing a Newton-type output least squares 
algorithm that reconstructs the admittivity distribution and the object shape simultaneously. The 
method is built in the framework of the complete electrode model and it is based on the Frechet 
derivative of the corresponding current-to- volt age map with respect to the object boundary shape. 
The functionality of the technique is demonstrated via numerical experiments with simulated mea- 
surement data. 
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1. Introduction. Electrical impedance tomography (EIT) is a noninvasive imag- 
ing technique which has applications, e.g., in medical imaging, process tomography, 
and nondestructive testing of materials [3j [31]. The objective of EIT is to recon- 
struct the admittivity distribution inside a physical body from boundary measure- 
ments of current and voltage. The most accurate model for EIT is the complete elec- 
trode model (CEM), which takes into account electrode shapes and contact impedances 
at electrode-object interfaces [6]. 

A real-life measurement setting of EIT typically contains more unknowns than the 
mere admittivity distribution: The exact electrode locations, the contact impedances 
and the shape of the imaged object are not necessarily known accurately. (As an exam- 
ple, consider a medical application where the body shape and the contact impedances 
vary from patient to patient.) These kinds of inaccuracies comprise a considerable 
difficulty for establishing EIT as a practical imaging modality since it is well known 
that even slight mismodelling can quite easily ruin the reconstruction of the admit- 
tivity [2j 0J |21] . The problems resulting from the aforementioned model uncertainties 
have partly been resolved in earlier works: Two alternative ways to handle unknown 
contact impedances have been introduced in [24, 33] , and fine-tuning the information 
on electrode positions has been considered in [8] . A brief review of the approaches to 
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tackling the problem with an unknown object boundary shape is given in the following; 
for a more extensive discussion, see [26] . 

Undoubtedly the most common way to treat problems resulting from an inaccu- 
rately known boundary shape is the use of difference imaging, where the alteration in 
the admittivity distribution is reconstructed on the basis of the difference between EIT 
measurements corresponding to two time instants (or frequencies) pQ. The method 
is based on the idea that the modeling errors are partly removed when difference 
data are used — given that the boundary shape remains unchanged between the two 
measurements. However, the difference imaging approach is highly approximative, 
because it relies on a linearization of the highly nonlinear forward model of EIT. 
Moreover, even if difference data are available, the boundary shape may also have 
changed between the measurements. This is the case, e.g., when imaging a human 
chest during a breathing cycle. A successful approach to coping with an unknown 
object boundary in absolute EIT imaging was suggested by Kolehmainen, Lassas and 
Ola [T9]l20]. Their method is based on allowing slightly anisotropic conductivities and 
on the use of sophisticated mathematical instruments such as quasiconformal maps 
and Teichmuller spaces. In [25] the so-called approximation error approach [18 was 
adapted to the compensation for errors resulting from an inaccurately known bound- 
ary shape in the framework of EIT. The approximation error method is based on 
the Bayesian inversion paradigm; the governing idea is to represent the error due to 
inaccurate modeling of the target as an auxiliary noise process. The (second order) 
statistics of the modeling error are approximated via simulations based on prior prob- 
ability models for the admittivity and the boundary shape. The application of EIT to 
imaging of human thorax was considered in [25] , where the approximated statistics of 
the modeling error were computed based on an atlas of anatomical CT chest images. 
In [26] , the method was further developed to allow the reconstruction of the boundary 
shape. See also [29] [30], where an optimization based technique was applied to the 
estimation of partially unknown boundary shape in process tomography applications. 

This work introduces an iterative Newton-type output least squares algorithm 
that tolerates uncertainties in the geometry of the imaged object. To be more precise, 
our aim is to include the estimation of the shape of the object boundary as a part 
of the reconstruction method. The required Frechet derivative of the measurement 
map of the CEM with respect to the exterior boundary shape is obtained with the 
help of domain derivative techniques stemming from [22] [12] [131 US 5 see a l so [IH] for a 
general theory of shape differentiation. However, unlike in [22 , 12 , 13 , 16 , the elliptic 
boundary value problem defining the derivative falls outside the standard i7 1 (0)- 
based variational theory due to Dirac delta type boundary conditions on the edges of 
the electrodes. This difficulty is tackled following the guidelines in [8 , where Frechet 
derivatives with respect to electrode shapes were considered, resulting in a well-posed 
'derivative problem' that is uniquely solvable in i^ 1_e (^), e > 0. 

Our approach is made computationally more tractable by introducing a dual 
method for sampling the i^ 1_e -regular shape derivative; in particular, it turns out that 
the reconstruction algorithm can be implemented without having to solve any forward 
problems with distributional boundary conditions. This observation is concretized by 
the numerical examples clearly demonstrating that the electrode measurements of 
EIT carry information on both the admittivity distribution and the object boundary 
shape. The numerical studies are based on simulated measurement data and carried 
out in three-dimensions, with the corresponding parameter choices founded on the 
Bayesian paradigm [T8] . 
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This text is organized as follows. Section [2] recalls the CEM and its fundamental 
properties. The main Frechet differentiability result is formulated in Section [3] and its 
proof is given in Section [4j Section [5] introduces the reconstruction algorithm, which is 
then tested numerically in Section [6j Finally, Section [7] lists the concluding remarks. 

2. Complete electrode model. Let Q C M n , n = 2 or 3, be a bounded domain 
and assume that its boundary dQ is an orientable C°°-manifold. We denote by a : Q — >• 
C nxn the electrical admittivity distribution of Q and assume that it satisfies the 
following, physically reasonable conditions [3]: 

a = a T , ReK-O^^I 2 , K ' ?l < C 2 \£\ 2 (2.1) 

for some constants Ci, C2 > and for all £ G C n almost everywhere in fi. 

Assume that the boundary dQ is partially covered with M G N \ {1} well- 
separated, open, bounded and connected electrodes {£Vn}m=i> i- e -> 

E m cdn, rn = 1 ..... M, and EjnE k = 0, 3 + (2.2) 

The electrodes are modelled as ideal conductors. The union of the electrodes is de- 
noted by E 1 = U m E m , and the frequency domain representations of the time-harmonic 
electrode current and potential patterns by the vectors / = [I m ]^ =1 and U = \Um]m=i 
of C M , respectively, where J m , U m G C correspond to the measurements on the rath 
electrode. Take note that the current vector /, actually, belongs to the subspace 

r 1 M 1 

Cf := [ Cl ,...,c M ]GC M Vc m = (2.3) 



due to to the current conservation law. The contact impedances (cf. [6]) that charac- 
terize the thin and highly resistive layers at the electrode-object interfaces are mod- 
elled by z G C M that is assumed to satisfy 

Re^>0, j = l,...,M. (2.4) 

According to the CEM EH], the pair (u,U) G ft 1 (ft) := (H 1 ^) © C M )/C, 
composed of the electromagnetic potential within Q and those on the electrodes, is 
the unique solution of the elliptic boundary value problem 

V • crVu = in ft, 

v-crVu = on <9ft \ E, 

u + z rn v • cfVu — Um on E mi ra = 1, . . . , M, 

^ • aVw dS = I mi ra = 1 , . . . , M, 



(2.5) 



for a given net electrode current pattern / G and with v = v(x) denoting the 
exterior unit normal of dft. The definition of ^(ft) as a quotient space emphasizes 
the freedom in the choice of the ground level of potential; in other words, one can 
never measure absolute potentials, only potential differences. 



The weak formulation of the CEM forward problem (2.5 ) is to find (u, U) G 1-1} (£1) 
that satisfies 28 



M 



B{(u, U), (v, V)} = trnVm for all (v, V) G ^(ft), (2.6) 
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where the sesquilinear form B: - H 1 (f2) x - H 1 (f2) — > C is defined by 

r M 1 r - 

B{(u,U),(v,V)} = / aVu- Vvdx+ V — / (U m - u)(V m - v) dS. (2.7) 
Jn m=1 z ™ J Em 

The form B is concordant with the natural quotient topology of H 1 (f2) (cf. [13 Corol- 
lary 2.6]), i.e., for all (u, U), (v, V) € U x {0) 

\B{(u, U), (v, V)}\ < dUK U)\\ nHQ) \\(v, V)\\ wm , 
ReB{(v,V),(v,V)} > C 2 \\(v,V)\\ 2 Him , 



where 



M 1/2 



\\(v, V)\\ ma) := mf {\\v - cf H1(Q) + \V m - c\ 2 } 



The unique solvability of (2.5) follows by combining the above estimates and the 
obvious boundedness of the antilinear functional on the right-hand side of (2.6) with 



the Lax-Milgram lemma [T5j |28] . This procedure also provides the estimate 

\\(U,u)\\ wm < C\I\, (2.8) 

where the constant of continuity C = C(Q,<j,z) can be chosen independently of the 
electrodes if it is assumed that 

min \E m \ > c (2.9) 

l<m<M 

for some constant c > (cf., e.g., [11, (2.4)]). In the rest of this work, we make the 
assumption (2.9) on the considered electrode configurations implicitly. 

An ideal measurement corresponding to the CEM provides the electrode voltages 
U G C M /C for some applied current pattern / G C^. For a given measurement 
setting {^,^,cr, z}, we thus define the measurement operator R: — )► C M /C by 

R.I^U. (2.10) 



Obviously, R is linear and bounded (cf. ( |2.5[ ) and ( |2.8[ )), with a constant of continuity 
that can be chosen independently of the electrode configuration under the assump- 
tion ( pH) ). 

To conclude this section, we note that for smooth a the interior potential has 
more regularity, namely 

v ■ (jVu\ e G H\E), v • a\/u\ dn G Hl /2 ~ e (dQ), u G H 2 ~ e (n)/C 

for all e > 0, as reasoned in [SJ Remark 1]. When appropriate, we emphasize the last 
statement by writing (u, U) g M 2_€ (1]) := (H 2 - e (Q) © C M )/C. 

3. Shape derivative. In this section, we introduce the derivative of the CEM 
measurement map with respect to perturbations of the object boundary dft. We begin 
by specifying how exactly the boundary is perturbed. 

For h G C^dft, R n ) we define 

F[h](x) = x + h(x), x G dQ, 
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and use the abbreviation dflh for the perturbed boundary, that is, 

dQ h = F[h](dQ) = {y eR n \ y = F[h](x) for some x e dft} . 

The open, origin-centered ball of radius d > in the topology of R n ) is denoted 

by B d , i.e., 

B d = {he C^d^R") | \\h\\ CHdW < d} . 

Following [10] , we introduce a special family of diffeomorphisms of R n to itself: 

Jo 1 = {F: R n -^R n | F-id £ C£(R n ,R n ) and F' 1 £ C 1 (R n ,R n )} , 

where Co(M n ,lR n ) denotes the space of continuously different iable vector fields that 
together with their partial derivatives vanish at infinity. In particular, when equipped 
with the natural norm, Co(lR n ,lR n ) is a Banach space (cf. [lOj p. 68]). The following 
proposition lists some fundamental properties of F[h] = id + ft for ft £ Bd with small 
enough d > 0. 

Proposition 3.1. There exists d = d(Q) > such that the following hold: 

(a) For every h £ Bd, ddh is the boundary of a bounded C 1 -domain Qh, and the 
mapping F[h] is a C 1 -diffeomorphism from dfl onto dQh'i 

(b) There exists an extension operator £ : Bd — > Co(IR n ,lR n ) such that 

£[h]\dn = K ||£[ft]||ci(R»,R»)< c ( n )\\ h \\ci(dn,R"<) 
and the extended mapping 

F[S[h]] = id + 5 [ft] 

belongs to for all h £ Bd- 

Proof The first part of the claim follows from an application of the implicit 
function theorem in local coordinates on dfl. The second part can be deduced, e.g., 
by first forcing h to zero in a tubular neighborhood of dQ and then using similar 
arguments as on page 78 of [10]. □ 

If there is no danger of a confusion, we abuse the notation by denoting the ex- 
tensions £[h] and F[£ [ft]] by the original symbols ft and F[h], respectively. Moreover, 
we assume implicitly that d > is as introduced in Proposition |3.1| 

Obviously, the measurement operator of the CEM may be considered as a map 
from B d x Cf to C M /C, i.e., 

R:(h,I)*->U[h], (3.1) 



where (u[h],U[h\) £ H 1 (^/ l ) is the unique solution of (2.5) when Q is replaced by 
Qh and the electrodes E m by := F[h](E m ) C <90^, m = 1,...,M. To make 
this definition unambiguous and to simplify the analysis that follows, we assume that 
a £ C 00 (R n ,C nxn ) with the bounds §2A§ satisfied everywhere in R n , i.e., that the 
admittivity distribution is defined in everywhere in R n — or at least in some proper 
neighborhood of Q. As a further simplification, we also assume that (in the three- 
dimensional case) the electrode boundaries dE m , m = 1, . . . , M, are smooth curves. 

We denote by ft r and h v the tangential (vector) and normal (scalar) components 
of ft £ Bd, respectively, that is, we have the (unique) decomposition ft = ft r + h v v. 
One might expect that it is enough to consider perturbations that belong to the 
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normal bundle of the boundary, i.e., ones that have vanishing tangential components. 
However, this turns out to be a false intuition, because tangential vector fields typically 
affect the measurement map in the 'first order' by moving the electrodes (cf. [8]) - 
even though they only define 'second order' perturbations of the object boundary dfl 
itself. 

Theorem 3.2. Under the above assumptions, the operator R: Bd x -» C M /C 
is Frechet differentiable at the origin with respect to the first variable, i.e., there exists 
a bounded bilinear operator R' : C 1 (<9^,R n ) xC M 4 C M /C such that 

lim l ] -— ^ \\R[h] - R[0] - R'h\\ c{C M^ M/c) = 

in C^dQ.R" 1 ). 

In the following, we will prove Theorem 3.2 in three dimensions, i.e. for n = 3, 
which is the more challenging case. The two-dimensional counterpart can be obtained 
by following a similar line of reasoning. 

The derivative R' of Theorem 3.2 can, in fact, be given explicitly. To this end, 
let H G C°°(dQ) be the mean curvature function defined so that it is positive if the 
surface turns away from the exterior unit normal, and consider the bounded surface 
divergence operator (cf., e.g., [7]) 

Div: [H s (dtt)]™ -> H s -\dn), s G R, (3.2) 

with the weak definition 

(Divv,ip) dn = -(v,Gmdip) dn , ip G C°°(<9ft), 

where Grad denotes the surface gradient (cf., e.g., [10 ]). We also introduce a family 
of distributions {S m }^ =1 C H-^ 2 - e (dQ), e > 0, defined through 

(S m ,v) dn = [ vols, v e H^^dQ), 

JdE m 

for m = 1, . . . , M. Notice that any v G HW+^dSi), e > 0, has a well defined 
restriction v\qe £ H e (dE) due to the trace theorem, and thus the definition of the 
family {S m }^ =1 is unambiguous. Moreover, we denote the characteristic function of 
E m C dQ by Xm-, m — 1 5 • • • •> M, and the unit exterior normal of dE in the tangent 
bundle of <9ft by vqe- 

With these tools in hand, let us consider the boundary value problem 

V • aW = in ft, 

M , m 1 



,M. 



v ■ aVu' - 2^ — (V - u')xm = fi+Z^ — (fcXm + fo5 m ) on dfl, 

m=l Zm m=l Zm 

I {U' m -u')dS = - ( f 2 dS- f hds, m = l, 

J Em J E rn JdEm 

(3.3) 

Here, the inputs h G H~ 1 / 2 ~ e {dQ), f 2 G ^^(E) and f 3 G H 1 ~ e (8E) are defined 
with the help of (u, U) G H 2_e (ft), i.e., the unperturbed solution of (2.5): 

fx = Div(h u (aVu\ dn ) T ), 
h\s m = K[(n - l)H(U m 
h\dE m = (h ■ v dB ){U m ~ u)\ 9 E m - 
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Notice that the claimed regularity of /i, fi and f% follows from (3.2) and (consecutive) 



applications of the trace theorem. It turns out that problem (3.3) is uniquely solvable 
in / H 1_e (^), and that the corresponding solution defines the Frechet derivative of 
Theorem IO 

Theorem 3.3. Under the assumptions of Theorem \3.2\ the boundary value 
problem (|3.3[) has a unique solution (u'[h),U'[h)) G H 1_e (£l), e > 0, for any 



h G C 1 (dQ,,W ri ). Moreover, the Frechet derivative of Theorem 3.2, i.e., R' 
C^dn.R 71 ) x Cf C M /C, is given by 

R' : (ft,/) i-> U'[h\. 



At first sight it may seem that Theorem |3.3| is not very practical as it defines the 
Frechet derivative of R with the help of a boundary value problem that falls outside 
the i^-based variational theory. Fortunately, there also exists a dual approach for 
sampling the shape derivative. 

Corollary 3.4. Let (u, U) G U 2 ~ e {yt), e > 0, be the solution of ([2^5} for some 
electrode current pattern IeCf. Then, for any (ft,/) G C^d^HT) x it holds 
that 

M r 

V(E>,J))J m = - h u (aVu) T -(Vu) T dS (3.4) 
„ i Jon 



i do, 

M 



du 



J2 — [ hJ(n-l)(U m -u)H 

M 1 f 

^ / ( h ' U dE)(U m ~ u)(U m - U) 



(Um-u) dS 



where (u,U) G H e (ft) is the solution of (2.5). 



4. Proof of the main result. Before moving on to prove Theorems 3.2 and 



3.3 and Corollary |3.4[ we give a brief summary of the variational technique on which 
the proof is based. A more complete reasoning in a slightly different framework can 
be found in [8j Section 5.1]. 

Let F = F[h] = id + ft G J-q, with ft G Bd, be as in the previous section. We 
introduce a pullback operator F* : i/ 1 (^/ l ) -» H 1 ^) defined by F*v = v o F\q; it 
is easy to see that F* is a linear isomorphism. A simple change of variables applied 



to the variational equation defining (u[h],U[h\) G H 1 (^/ l ), cf. (3.1), shows that the 
difference of the pullback pair (F*u[h], U[h]) G H 1 ^) and the unperturbed solution 
(u,U) = (u[0],U[0\) G H 1 ^) satisfies (cf., e.g., [16]) 

B{(F*u[h] - u, U[h] - U), (v, V)} 

= [ (p-<r*[h]) V(F*u[h]) • Vvdx 

M _ 

S— / (U m [h}-F*u[h})(V m -v)(l-\JncF\)dS (4.1) 



m=l 



for all (v,V) G H 1 (^). Here, the pullback admittivity <r*[ft] is defined as 

a*[h] = \J F \J-\F*a)(J- 1 ) T , 
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Jf is the Jacobian matrix of F, \Jf\ is the absolute value of its determinant, and Jac F 
is the surface Jacobian determinant of the restriction F\qq : dft — » dQh- Moreover, 
it follows from the perturbation analysis in [T2l [13] that modulo 0(||^||^i( R n r it 
holds that 

a - a*[h] = o-Jl + J h a - (ft • V + V • ft)a, (4.2) 

l-|JacF|= -(n-i)Hh„-Divh T , (4.3) 

where ft • Va is defined as the matrix (ft • Vcr^)^ J=1 . In consequence, in order to 
estimate the difference (F*u[h] — u,U[h] — U), it seems reasonable to consider the 
bounded sesquilinear functional A : C^R^R") x U 1 ^) C defined by (cf. 
(19)]) 

A[h](v,V)= [ (aJl + J h o--(h-V + V -h)(j)Vu-Vvdx 
Jn 

M i r - 

— / (U m -u)(V m -v)((n-l)HK + Dwh T )dS, 

m= l Z ™ JE m 

and the corresponding ft-parametrized variational problem 

B{(w[h],W[h]), (y, V)} = A[ft](v, V) for all (v, V) G W^fi), (4.4) 

which has a unique solution in due to the Lax-Milgram lemma. The following 

proposition is a straightforward variation of [8, Proposition 5.4]. 

Proposition 4.1. Let (n, U) G H 1 ^) and (n[ft], 17 [ft]) G ^(flh) be as defined in 
Section\^and (w[h],W[h]) G H 1 ^) tfie unique solution of ( |4.4| ). Tften, £fte estimate 



\\(F*u[h] - u, U[h] -U)- (w[h],W[h])\\ n i (n) < C\I\\\h\\l 1{my 



holds with a constant C > £fta£ can 6e chosen independently of I G and ft G S^. 

Although the map ft i->> W[ft] is a first order approximation of ft ^ (R[h] — 
R[0])I = U[h] — U[0] around the origin, it does not provide a satisfactory definition 
for the Frechet derivative ft \-> R'[h\. Indeed, although ft \-t W[h] is clearly linear with 
respect to the extension ft = £ [ft] G Cq (W 1 , R n ) , it is not self-evident that the s ame 



also holds for the original perturbation ft G C 1 (dft,lR n ) as required by Theorem 3.2 
Moreover, from the computational view point, the extension of ft to the whole of R n 
is a nuisance that one wants to avoid. 

To get rid of this problem, we proceed as in [8] and modify the first component of 
(it; [ft], W [ft]) in an appropriate way. This procedure involves including a directional 
derivative of the interior potential component of the unperturbed solution (u,U) G 
l~L 2 ~ e (yt) as an argument of the sesquilinear form B: / H 1 (f7) x 1-1} (Q) — » C. Since the 
derivatives of u are merely in i7 1_e (ft) such analysis cannot be carried out without any 
modifications. For this reason, we introduce a sequence of smooth approximations for 
u G i7 2_e (ft)/C, and subsequently also for (w[h], W[h]) G H 1 ^). As in [8j Subsection 
5.4], we may pick a sequence (u^\U^) G (C°°(Q) © C M )/C, j = 1, 2, . . . , such that 

V-crW^=0 in ft and lim (u u \ U U) ) = (n, U) in H 2 " e (ft). (4.5) 

Moreover, we define (w^[h],W ij) [h]) G V}(Vt) to be the unique element of V}(Vt) 
that solves the variational problem 

B{(w^[h],W ij) [h]),(v,V)} = Aj[h](v,V) for all (v,V) G H 1 (ft), (4.6) 
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where the antilinear functional Aj[ft]: — » C is defined via replacing (u,U) by 

(u^\U^) in the definition of A [ft]. Through a slight variation of the argument in 
the proof of [SJ Lemma 5.7], one easily obtains that 

lim (w&[h], W {j) [h}) = (w[h],W[h]) in H\n). (4.7) 

3 -too 

We proceed by defining the 'augmented interior derivatives' by 
w[h]=w[h]-h-Vu, w U) [h]=w U) [h]-h-Vu u \ j = 1,2, ... . (4.8) 
Due to ( [45] ) and ( pj~7| ), it follows that 

lim (wW[h],wW[h]) = {w[h],W[h]) in H 1 ' 6 ^) (4.9) 

j->>00 

for any e > (cf. [23 ). In the following, we will show that (w[h], W[h]) G H 1_e (^) is 



the unique solution of (3.3) for ft G Bd- In particular, the pair (w[h], W[h}) turns out 
to be independent of the extension of ft G C 1 (<9^,R n ) to the whole of R n . 

Proof of Theore ms\3.£\ and \3.3\ In the first part of the proof, we show that the 
derivative problem ( |3.3[ ) is uniquely solvable, with the corresponding solution being 
the above constructed pair (it) [ft], W[ft]) G H 1_e (^) if ft G Bd- First of all, we note 
that the uniqueness of the solution can be proved in exactly the same way as in the 
first part of the proof of [HI Theorem 6.1], which means that we can focus solely on 
the existence. In case that v • ft = and ft G Bd, the fact that (?2)[ft], W[ft]) is a 



solution to (3.3) follows from the second and third parts of the proof of [8[ Theorem 



6.1]. Due to the linearity of the right-hand side of (3.3) with respect to ft, the unique 



solvability of (3.3) thus follows if we are able to show that (w[h], W[ft]) is a solution 
also if ft = h v v G Bd- (For a general ft ^ Bd the unique solution of ( |3.3[ ) can then 
be obtained by rescaling the corresponding solution for h/c G Bd with large enough 
c>0.) 

(1) Assume that ft = h v v G Bd and let us consider what kinds of variational 
equations (w^[h],W^[h]) G V 1 ^) satisfies. For now, let ip G C°°(n) and V G C M 
be arbitrary. Recalling first and the definition of (w^[h},W^[h\) G U l (Q), 

and then using standard vector calculus, the first part of ( |4.5[ ) and the divergence 
theorem, we obtain (cf., e.g., [T6] ) 

} dTp 

on 

M 

z 

m=l 

M 



B{(w^[h],W^[h]),{(p,V)} = h v v (aVu^^ - (cjVu {j) -Vlp)v)dS 

Jf)0 \ Oh> J 

A l r du& - 

V — / K^-(V m -lp)dS (4.10) 

m=l Z m J Em 
M r 

^ L 

71=1 Zm JE ' 



KH(U$ -u^){V m -Tp)dS. 



By dividing Vu^ and V(f into tangential and normal components on dfl, the first 
integrand further simplifies as 

h v v • (aVu^^ - (aW j) • VTp)v) = -h v (o-Vu^ j) ) T • (Vp) r . 
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Due to (4.5) and the regularity of the unperturbed solution (u,U) G V 2 e (^), we 
may take the limit j — >• oo in (4.10), yielding 



lim B 



{(w^[h],W^[h]),(cp,V)}=- [ h v {crVu) T .{Vy) T dS 

Jon 

M i r 8 

+ E— / h^{V m -Tp)dS (4.11) 

m=l Zm Je ™ OV 
M -If 

- / h * H ( u ™ - u){V m - Tp) dS. 



To prove the first equality of (3.3), let V = and ip G Cq°(Q) be arbitrary. 
According to the definition of the sesquilinear form £?, the identity ( |4.10 ) and the 
definition of distributional differentiation (cf., e.g., [9]), it holds that 



(V-<rVw W) [ft],<p) n =0. 



(4.12) 



As the elliptic differential operator V • <rV: H 1 ~ e (Q)/C — »• i^ _1_e (£7) is continuous 
for any e G R such that e — 1/2 ^ Z [23[ Chapter 1, Proposition 12.1], we may pass 
the limit inside the brackets of (4.12). Consequently, V • crVw[h] = is satisfied in 
the sense of distributions in Vt. 

Assume next that cp G C°°(Q) and let still V = 0. The (generalized) Green's 
formula (cf. p. 382, Corollary 1]) and fl4 12| ) indicate 

B{(wW [ft] , [ft] ) , (<p, 0) } = (v • aVw® [ft] , ^) afi 

M 

m=l ^ m 

Moreover, according to [23j Chapter 2, Theorem 7.3], the Neumann trace map v \-> 
v • <j\7v\qq is well-defined and bounded from the closed subspace 

{v G fr 1_e (fi)/C | V • aVv = 0} C H^tfy/C 

to i^ _1 / 2_e (9^). Thus, ( |4.9[ ) and the trace theorem give 

.lim (</>,())} = (iz-trVtiiW,^ 

M 



J E. 



On the other hand, by (4.11) it also holds that 

t{(wW[h],wW[h]),{<p,0)}=- [ h v {aVu) T -{VTp) T dS 

Jon 

M 



lim B\ 



m=l 
M 

E 



J E, 

n — 1 



hv^—pdS 



h v H(U m — u)TpdS. 
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As (V^) r = Grad^ on dQ (cf., e.g., [7 ) and C°°(n)\ dQ is dense in H a (dQ) for 
any sGl, this proves that (w[ft], W[ft]) satisfies the second equation of (3.3), since 
h • VdE — by assumption. 

To prove the remaining, third condition of (3.3), let V be the rath coordinate 
vector and choose tp = 0. Using the definition of £?, (4.9) and (4.10), we conclude 
that 



(W m [h]-w[h])dS- 



l)H(U m - u) 



du 
dv 



(4.13) 



Since ra was chosen arbitrarily, it follows that (w[h], W[h]) G V 1 e (0) is a solution 
to Q for ft G Bd- 

(2) Let us then prove that the mapping 



R' : (ft,/) ^ U'[h], C^dn.R" 1 ) x Cf -> C 



M , 



really defines the Frechet derivative of Theorem |3.2| as claimed in Theorem 3.3| First 
of all, it is easy to see that R' is bilinear since the right-hand side of ( |3.3 ) depends 
bilinearly on ft G C 1 (9D,lR n ) and (u,U), and the unperturbed solution (u,U) itself 
depends linearly on the applied current pattern. Moreover, due to Proposition |4.1| 
and since U'[h] — W[h] for ft G Bd by the first part of the proof, we may estimate as 
follows: 

\\U[h] -U-U'[h]\\ C M /c < C\I\\\h\\ 2 CHdClt1tn) , h e B d , 

which completes the proof as C > can be chosen independently of / G Cf and 
heB d . □ 
We complete this section by providing a proof for Corollary |3.4| 
Proof of Corollary \3.4\ As in the previous proof, it is enough to consider small 
ft in the normal bundle of dfl, i.e., ft = h v v G Bd, by the virtue of the linearity 
of the claimed sampling formula with respect to ft and the fact that for tangential 
perturbations the assertion follows through the same line of reasoning as [8j Corollary 
4.2]. 

Let I G Cf and U) G be as in Corollary 3.4 and consider ft = h v v G Bd- 

Due to Theorem [3^1 the fact that (V [ft], U'[h) = (w[h],W[H\), the limit {49| and the 
variational formulation (|2.6|) corresponding to the current pattern /, it holds that 



M 



M 



V/ m (i?'(M))m= lim yZlmW^[h} = lim B{(w^[h],W^[h]),(u, U)}. 



On the other hand, following the same line of reasoning as in (4.11) — and approxi- 
mating u by a sequence of smooth functions {(fj} — we obtain that 



lim B{(wW[h], W U) [h]),(u, U)} 



M 



h u (aVu) r ■ (Vu) r dS 



m=l 
M 



u) dS 



V ^-^ / KH{U m - u)(U m - u) dS, 
m=1 z ™ J Em 



which is the normal bundle version of (3.4) and thus completes the proof. 



□ 
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5. Algorithmic implementation. In this section, we introduce our numeri- 
cal algorithm for the simultaneous reconstruction of the admittivity distribution and 
the object boundary. It is assumed that the object of interest Q C M 3 is a cylin- 
der D x (0, ho), where D C M? is a simply connected and bounded cross-section 
shape, and ho > the known height of the body. The electrodes are of the form 
E m = 7m x (0, ho), with each j m being a connected part of dD with a known length, 
i.e., the electrodes are rectangular, homogeneous in the vertical direction and assumed 
to be of a known width and the same height as the object itself. In particular, if the 
admittivity distribution were also homogeneous in the vertical direction — as it is 
in our numerical experiments — , the measurement setting could be modelled by a 
two-dimensional forward problem. Be that as it may, we carry out all numerical com- 
putations in three dimensions in order to demonstrate the feasibility of our method 
in a realistic framework. Moreover, we only consider real-valued and isotropic elec- 
trical admit tivities, i.e., a : Q — >• R + . Note that the above assumptions on the target 
are made only for the sake of simplicity; the generalization of the algorithm to more 
general three-dimensional settings is conceptually straightforward. 

In the following three sections we outline the ideas behind our reconstruction 
method, but do not discuss all details about, e.g., the form of the smoothness prior 
for the admittivity; see, e.g., [TTJ [18] for more information. 

5.1. Parametrization of the unknowns. In many practical situations the 
examined body has a star-shaped cross-section. In consequence, we search for the 
unknown boundary dD as a C°° -curve parametrized by 



la (4>) 



N 

COS ( 



ao + ^^( a j cos j(J) + <x/+at sin j</>) 

3 = 1 



sine 



a , • • -,OL2N € ^, (5.1) 



where <\> is the polar angle and the coefficients a = [ao,..., «27v] T £ 
sumed to be such that the curve does not intersect itself. Let D a denote the bounded 
set of R 2 with dD a = 7 a ([0, 27r]) and furthermore define Q a = D a x (0, ho). As 
it is assumed that the width of the (rectangular) electrodes is known, we may thus 
parametrize them by their initial polar angles m , m = 1, . . . , M, in the counterclock- 
wise direction. The vector containing these angles is denoted by = [0\, . . . ,0m] T - 
We assume that the electrodes are numbered in the natural order, that is, the terminal 
angle of an electrode precedes the initial angle of the following one. 



Approximate forward solutions to (2.5) in Q a are computed by a finite element 



method (FEM). The FEM solver used in this work is an adaptation of the implemen- 
tation in [32 . In the FE scheme, we discretize the computational domain Q a into 
tetrahedrons and approximate the distributions of admittivity and potential in piece- 
wise linear and quadratic bases, respectively. In our reconstruction algorithm, the 
geometric parameters a and change iteratively and consequently Q a and its FEM 
mesh also change at each step. In order to fix the admittivity discretization indepen- 
dently of such deformations, we pick a sufficiently large cylinder E = B x (0, ho), with 
a discoidal base B inside which we let the cross-section D a evolve. Given this back- 
ground cylinder, we look for admittivities of the form ^ k a^Lp^, where {cr/e} C (0, oo) 
and {(ft} is the piecewise linear basis related to a fixed tetrahedral mesh of E. The 
admittivity values are transformed between the fixed 'reconstruction mesh' of E and 
the varying ones in Q a via linear interpolation. 

5.2. Bayesian framework. Although our reconstruction algorithm, which will 
be introduced in Section |5.3[ cannot be considered purely Bayesian, its underlying 
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motivation is statistical, and thus the basic ideas behind Bayesian inversion are out- 
lined in the following. In the Bayesian approach, all quantities are considered as 
random variables with some assumed prior probability distributions. Combining the 
information from the prior with the measurement data, one gets the updated posterior 
distribution for the parameters of interest [18] . 

Let {I^}^! 1 be a basis of . The voltages measured at the electrodes on the 
boundary of Q a are modelled as 

y(j) = u( j Xa,a,6) + r)( j) Gl M , j = l,...,M-l. (5.2) 

Here, (u^\a,a,9),U^\a,a,9)) G H 1 ^) © is the solution of in the do- 



main Q a with a real- valued admittivity <r, when the net electrode current pattern 
1^ is injected through the M electrodes parametrized by their initial polar angles 
6 G R M . Notice that we have fixed the ground level of potential by requiring that 
, 9) has vanishing mean. The components of the noise vector rf^ G R M are 
assumed to be independent realizations of zero mean Gaussian random variables. To 
simplify the notation, we pile the electrode currents, potentials and noise vectors into 
arrays of length M 2 — M, i.e., employ the shorthand notations 

x=[(/( 1 )) T ,...,(/^- 1 )) T ] T , 
v = [(W 1 )) T ,...,(W^- 1 )) T ] T 

r / =[(^...,(^- 1 ))Y [ ' } 

U(a, a, 9; X) = [U™ (a, a, #) T , . . . , U^ M ~^ (a, a, #) T ] T . 

This allows us to write the noisy measurement model as 

V= U(a,a,0;X) + r] eR M2 ~ M . (5.4) 

The discretized admittivity is given a homogeneous Gaussian smoothness prior 
with a covariance T a and a positive homogeneous mean a*; for more details about 
smoothness priors, see [18]. To include control over the geometric information, the 
coefficients a are provided with a Gaussian prior with a mean a* G M 2Ar+1 and a 
covariance matrix r a = diag(ag, . . . , cl^n)^ wnere 

fr s a, l = 0,...,N 

a* = < , x (5.5 
3 \(l-N)- s a, / = 7V + 1,...,27V. v 7 

By adjusting the parameters s, a > 0, one may tune the prior assumption on the 
regularity of the object boundary. The electrode initial polar angles are also given 
a Gaussian prior density with a mean 0* G R M and a diagonal covariance matrix 
Yq = t 2 I, where r > is the corresponding standard deviation. As a result, for a 



measurement V of the form (5.2), a maximum a posteriori (MAP) estimate is obtained 



as a minimizer of the Tikhonov-type functional 

$(a, a, 9) := (W(<7, a, 6>; X) - VfT' 1 (U(a, a, 6>; X) - V) + (a - a*) 1 ^- 1 (a - a*) 

+ (a - a*) T T-\a - a*) + r~ 2 \9 - 9* | 2 , (5.6) 



where is the diagonal noise covariance matrix; see [T8] . 
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5.3. The (quasi-Bayesian) iterative algorithm. According to our experi- 
ence, the EIT measurements modelled by the CEM are typically more sensitive to 
the exterior boundary shape and electrode locations than to the internal admittivity 
distribution. As a consequence, it seems to be computationally advantageous to first 
fix a crude constant approximation for the admittivity distribution, then use a (de- 
terministic) iterative scheme to come up with a relatively good model for the object 
boundary and electrode positions, and finally use these preliminary estimates as the 



prior expectations in the to-be- minimized MAP functional (5.6). It should be empha- 
sized that such an initialization of the means makes our algorithm strictly speaking 
non-Bayesian, since the choice of the priors should be independent of the data. Be 
that as it may, according to our experience, such a preliminary step leads to faster 
and more reliable convergence. (We do not claim, however, that this kind of two-step 
implementation is the only feasible choice.) 

To be more precise, we first choose the covariance matrices T^, r a , T a and Tq 
according to the assumed prior information on the variation of the corresponding 
parameters, pick an initial guess (a^°\0^) for the measurement geometry (corre- 
sponding to some disk-shaped cross-section in all of our numerical studies), and fix cr* 
to be the constant admittivity that minimizes the output least squares part of <£, i.e., 

(W(<7, a, 0; X) - VfT' 1 (U(a, a, 0; X) - V) 

when (ce, 0) = (a^°\0^). Subsequently, the following two-stage scheme is employed: 

First stage of the algorithm: choosing the prior means. We apply a 
Levenberg-Marquardt type method in order to choose the geometry parameters that 



are used as the prior means (a*, 0*) when of (5.6) is minimized simultaneously with 
respect to all of its variables in the second stage of the algorithm: 

1. Fix a — o~* in ( |5.6[ ), consider $ as a function of only two variables a and 0, 
and set (a*,0*) = (a(°),0(°)). 



2. Calculate the Gauss-Newton minimization direction for $(a, 0) of (5.6) at 
(a,0) = (a*,0*); see, e.g., [27]. 

3. Minimize $(a, 0) over the line passing through (a*, 0*) in the Gauss-Newton 
direction by the Golden section line search. Redefine (a*,0*) to be the ob- 
tained minimizer. 

4. If satisfactory convergence is achieved, terminate the iteration. Otherwise, 
return to step 2. 

This part of the reconstruction algorithm is stable and does not, in particular, seem 



very sensitive with respect to the choice of the covariance matrices in (5.6). 



Second stage of the algorithm: finding the MAP estimate. After the 
prior means (a*,a*,0*) have been chosen in the earlier parts of the algorithm, the 
final stage consists of minimizing <1> of ( |5.6[ ) by the Gauss-Newton algorithm: 

1. Set k = 1 and (^ k \a^\0^) = (cr*,a*,0*). 

2. Calculate the Gauss-Newton direction for 3>(cr, a, 6) of ( |5.6| ) at (cr, a, 0) — 
(ffO.aO.jW). 

3. Minimize $ over the line passing through (<jW,aW,0W) in the 
Gauss-Newton direction by the Golden section line search and define 
(cr^ +1 \ a( k+1 \ 0( /c+1 )) to be the obtained minimizer. 

4. Unless satisfactory convergence is achieved, increase k by one and return to 
step 2. 
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For the computation of the needed Gauss-Newton directions, one needs the Ja- 
cobian of U(<j,a,0;I) with respect to <r, a and 0. By the Jacobian with respect to 
a we mean the one with respect to the coefficients of the piecewise linear basis in 
the 'reconstruction cylinder' E introduced in the last paragraph of Subsection 5.1. 
(Notice that the coefficients of the basis functions supported outside Q a do not play 



a role in U(a, a, #;Z), but they do affect the last term on the first line of (5.6).) For 



the estimation of the derivatives with respect to a and 0, we refer to [17] and [8], re- 



spectively. By the dual relation (3.4), the Jacobian with respect to a can be sampled 
via trivial linear algebra (a change of basis) after evaluating the expressions on the 
right-hand side of (|3.4|) for each triplet 



h(<f>) = ^ l \<p) 



COS ( 

sine 



(u, U) = (u^\ /7 (i) ), (fi, U) = (vV\uW) 



over the indices I = 0, . . . , 27V and i, j = 1,...,M- 1. Here, (u^\U^) = 
(uV\a,a,0),U( j \a,a,0)) is the solution of (ESJ) for I = and the setting 



parametrized by (a, a, 0), and ipi((/)) = coslcj) if I < N and ^(0) = sin(Z — AT)0 
when I > N + 1. We emphasize that one needs not solve any extra forward problems 
for this procedure since at each iteration step all the pairs (u^\a, a, 0), U^\a, a, #)), 
j = 1, . . . , M — 1, must be computed already for evaluating the functional $ of ( |5.6[ ). 
By the assumption that the electrode width is known, for any given electrode the 
terminal polar angle is a smooth function of the initial one and a. This functional 
dependence can be written explicitly by employing the arc length formula for the 



parametrization (5.1), and this information can then be included in the Jacobians 
with the help of the Leibniz rule and the chain rule for the total derivative. 

Remark 5.1. If one chooses to skip the first stage of the above introduced algo- 
rithm and use the (simple) initial guess (a^°\ 0^) as the prior mean for the geometric 
parameters in the second stage, with suitable parameter choices the reconstructions for 
the numerical experiments of the following section typically remain qualitatively the 
same, but the convergence slows down considerably. What is more, in practice the 
initial guess (a^°\0^) for the measurement geometry is often more accurate than 
the ones we employ in our numerical experiments, which further reduces the practical 
relevance of the first stage of the algorithm. 

6. Numerical experiments. Our main aim is not so much to compare the 
functionality of our method with reconstruction techniques presented elsewhere, but 
to make an 'internal' comparison between three cases: 

(i) the measurement geometry, i.e. the object shape and the electrode locations, 
is known; 

(ii) the measurement geometry is known inaccurately but this is not taken into 
account in the algorithm; 

(iii) the unknown boundary shape is estimated simultaneously with the admittiv- 
ity distribution. 

We will demonstrate that (i) and (iii) give comparable results, while the quality of 
reconstructions for (ii) is intolerably bad. 

We present three numerical experiments, in each of which M = 16 identical elec- 
trodes of known width are attached to the object of interest. We assume to know the 
contact impedances and choose the values z m = 1, m = 1, . . . , M. The first exper- 
iment, though a bit impractical, acts as an initial probe to test the functionality of 
the computed Frechet derivatives: we apply (the first part of) our algorithm to the 



16 



J. DARDE, N. HYVONEN, A. SEPPANEN AND S. STABOULIS 



shape estimation of a target object with a known homogeneous admittivity distribu- 
tion. In the second experiment we consider a simple shape (an ellipse) and a smooth 
admittivity distribution. In the last experiment the object shape is moderately com- 
plicated and the admittivity phantom consists of inclusions of constant admittivity 
in a homogeneous background. 

Let q be the to-be-reconstructed admittivity and suppose the pair (/?, $) provides 
a parametrization of the target measurement setting in the sense of Section [5TT] To be 
quite precise, the latter statement is a bit ambiguous because none of the considered 
target shapes dD can be given in the form ( |5.1[ ) with a finite TV, but for ease of 
notation we have decided to allow here an 'infinite' shape parameter vector (3. For 
each experiment we simulate the exact data U(s : /3, using the input current basis 

= Gl _ e . £ j = 2, . . . , M. where e^ is the jth Euclidean basis vector. Notice 
that there is no danger of an inverse crime because a new finite element mesh for the 
approximate domain Q a is generated at each iteration of the reconstruction algorithm, 
and these meshes differ considerably from the mesh of the target object = used 
for the data simulation. 



The actual noisy measurement realization V is formed via (5.4), with (a, a, 9) = 
(5, /?, by picking a particular noise component r]m\ m = 1, . . . , M, j = 1, . . . , M— 1, 
from a zero mean Gaussian distribution with the variance 



0.01 2 |[^)( ? ,/^)| 2 



0.001 2 max \uj?\<;, /?, 0) - U"> (?,/?, 0)| 2 
i<fc,;<M K 1 



Td)l 



(6.1) 



Here, the relation between U(<;, f3,$;T) and Um\s,f3,$) is as in (5.3). Sections 6.2 
and 6.3, where the cases (i-iii) are compared, work with fixed realizations of the noise 
vector rj in order to allow a fair comparison. 



For further justification of the noise model (6.1), see [15] , but anyway notice that 



(6.1 ) corresponds to more than one percent of relative noise in the absolute data, which 



is a substantial amount for an EIT problem. In each numerical experiment we assume 
to know the covariance of the measurement noise, i.e., we use the diagonal covariance 



matrix defined by the noise model (6.1) as in (5.6). We do not elaborate on the 
choice of Y a in (5.6) in further detail; it is built based on the proper (informative) 
smoothness prior proposed in [18 , reflecting the a priori assumption on the spatial 
variations of the admittivity. 



6.1. Known homogeneous target admittivity. Figure 6.1 shows the results 
obtained when the first stage of our algorithm is applied to reconstructing the bound- 
ary shape and electrode locations for an object with a known homogeneous admittivity 
distribution <; = 1. While this situation has minor practical relevance, it serves as a 
test of the computational techniques for obtaining the derivatives with respect to a 
and 6. The cylindrical target object is Q = D x (0,/io) with ho = 1, and the curve 
dD is parametrized by 



7W = 



3.3 



(2.2 2 cos 2 0+i.5 2 sin 2 0)1/2 



Lie 



-(0-tt) 6 



•0.88 cos^sin(-2(/>) 



COS ( 

sin c 



The width of the M = 16 identical electrodes on dD x (0, ho) is 0.3 and their initial 
polar angles are of the form r d rn = 2i\(m — 1)/M + e m , m = 1, . . . , M. where e m are 
independent realizations of a normally distributed random variable with zero mean 
and standard deviation 0.1. 

Since the admittivity is known a priori, we do not estimate a* as explained in the 



beginning of Section 5.3 , but fix it to be identically 1. The first stage of the algorithm 
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(a) The iterates. (b) Reconstructed geometry. 



Fig. 6.1: Retrieval of an unknown cross-section shape from noisy simulated data, 
(a) The five iterates with the final one plotted with solid line, (b) Comparison between 
the exact shape (red solid) and the retrieved one (blue dashed). 



in Section 5.3 is then run with N = 15 and M = 16; in this case the final value of 



(a*,0*) describes the reconstructed measurement geometry. In the construction of 



the prior covariance T a we use (5.5) with the selection a = s = 1 and set Tq = r 2 I 
with r = 27r/M, but the algorithm does not seem to be very sensitive with respect to 
these choices. Figure |6J^ a) shows the required five iteration steps, one of which is the 
eventual reconstruction, obtained by choosing the initial guesses = [2.7, 0, . . . , 0] T 
and 0$ = 27r(ra — 1)/M, m = 1, . . . , M — 1. The final iterate is drawn with solid 
line and the others with dashed line. In Figure |6J^ b), the target curve 7 (red solid) 
is compared with the retrieved one (blue dashed). 

The algorithm was run with several different target objects and in all cases the 
results were qualitatively similar to what is illustrated in Figure |6.1[ given that the 
examined shapes were not too complicated: If there were fine structures on a scale 
smaller than the electrode width, the results were poor. Further, the simpler the 
geometry, the faster the convergence was. It was also observed that the number 



of coefficients in (5.1) should not be too large, at most about N = 15. A high 
number of coefficients results in unstable reconstructions and absurd shapes, with the 
performance of the algorithm slowing down. 

6.2. Smooth target admittivity. In the second experiment, we apply the 
(whole) simultaneous reconstruction algorithm to data corresponding to a relatively 



simple target shape and a smooth admittivity distribution illustrated in Figure 6.2 ^a). 
The shape of the target object is Q = D x (0, /io), ^0 = 0.5, where D is an ellipse with 
major and minor semi- axes 2 and 1.5, respectively. The admittivity is homogeneous 
in the vertical direction, which allows us to only consider cross-sections in the visu- 
alizations. The electrode positions are chosen in the same manner as in the previous 
example, with the constant electrode width being such that two fifths of dD x (0, ho) 
is covered by the electrodes. 

In the reconstruction process we seek for a parameter triplet (a, a, 0) G M x x M N x 
R M with N = 7 and M = 16. Here, K is the number of nodes in the (fine enough) 
discretization of the background cylinder E = B x (0,/io), with B chosen to be an 
origin-centered disk of radius 3. As the initial guesses, we use a (°) = [2,0,...,0] T and 
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(a) Phantom. 



(b) Incorrect geometry. 




(c) Correct geometry. 



- " 1.6 

- " 1.4 

- " 1.2 

jo.8 
Ho.6 




(d) Simultaneous retrieval. 



Fig. 6.2: Experiment with a simple target shape and a smooth admittivity; each 
image represents a cross-section of the corresponding phantom/reconstruction that is 
(almost) homogeneous in the vertical direction, (a) Phantom used in data simulation, 
(b) Reconstruction corresponding to an incorrect fixed geometry, (c) Reconstruction 
corresponding to the exact geometry, (d) Simultaneously reconstructed admittivity 
and measurement geometry. 

6$ — 2ir(m — 1)/M, m = 1, . . . , M — 1. The prior covariances are constructed by 
selecting a = 0.1, s = 1 for r a of (5.5) and Tq = r 2 I with r = 2tt/M; see the fourth 



paragraph of Section [6] for an explanation about the choice of T a and T^. 

The reconstruction in Figure [6^ b) was obtained by applying the second stage of 
the algorithm in Section [573] with respect to a to the setting where the last two terms 
of (5.6) are deleted and (a, 6) is fixed to be the initial guess (a^°\ this approach 



corresponds to ignoring the incompleteness of the information on the measurement 
configuration and assuming stubbornly that the cross-section of the target object 
is a disk with uniformly distributed electrodes on its boundary. The reconstruction 



corresponding to the precise knowledge of the geometry is depicted in Figure 6.2 'c); it 
was obtained in the same manner as the one in Figure 6.2 b), except this time around 
the geometry parameters (a, 0) = (/?,#) were fixed at the values describing the target 
configuration used in the simulation of the measurement data. Figure 6.2 ^d) visualizes 
simultaneous retrieval of the admittivity distribution and the measurement geometry; 
this reconstruction c orre sponds to application of the whole two-stage reconstruction 
algorithm of Section 5.3 , starting from the initial guess (a^\0^) defined above. 
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From Figure 6.2 ^b), it is obvious that ignoring the incompleteness of the informa- 
tion on the measurement configuration results in severe artefacts in the admittivity 
reconstruction close to the object boundary. On the other hand, a comparison of Fig- 
ure [6]2^d) with Figures [cT^c) and |6.2[ b) demonstrates that the simultaneous retrieval 
of the admittivity distribution and the measurement setting provides a qualitatively 
similar reconstruction as knowing the exact geometry to begin with, and a far better 
one than altogether ignoring the inaccuracies in the geometric information. 

6.3. Piecewise constant target admittivity. In our last experiment, we 
consider the target object illustrated in Figure [6^ a). It is characterized by Q = 
D x (0, ho), ho = 0.5, with dD parametrized by 



7W = 



(1.5 2 < 



2 2 : 



+ 0.75e- ((/) - 7r)6 +O.6cos0sin(-2(/>) 



COS ( 

sin c 



The corresponding admittivity distribution, which is homogeneous in the vertical 
direction, consists of a homogeneous unit background and two embedded inclusions 




*4 







(a) Phantom. 



(b) Incorrect geometry. 




(c) Correct geometry. 



(d) Simultaneous retrieval. 



Fig. 6.3: Experiment with a complicated target shape and a piecewise constant 
admittivity; each image represents a cross-section of the corresponding phan- 
tom/reconstruction that is (almost) homogeneous in the vertical direction, (a) Phan- 
tom used in data simulation; the admittivity of the inclusions is 10. (b) Reconstruction 
corresponding to an incorrect fixed geometry, (c) Reconstruction corresponding to the 
exact geometry, (d) Simultaneously reconstructed admittivity and measurement ge- 
ometry. 
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with the constant admittivity level 10. The target electrodes are of equal width, they 
cover two fifths of dD and their locations are chosen as in the previous examples. 

pK . 



In this case we consider <1> of ( 5.6 ) as a function of (a, a, 6) G R K x R 7V x R M , with 
N = 15, M = 16 and K being the number of nodes in the mesh for the background 
cylinder E = B x (0, ho). Here, B is once again a disk of radius 3 centered at the 
origin. We assume the same prior information as in the previous example: T a is as 

r 2 I with r 



in (5.5) with a = 0.1, s = 1 and Yq 



the iterative reconstruction algorithm of Section 



3 (o) 



5.3 



2tt/M. The initial guesses for 
[1.5,0,...,0] T and 



are a 



(o) 



27r(ra - 1)/M, m 



,M -1. 



The results are illustrated in Figure 6.3 , with the subimages organized as in 



Figure 6.2 of the previous section. The reconstruction shown in Figure |6.3[ b) was 
obtained by ignoring the incompleteness of the information on the geometry, i.e., ap- 
plying the second stage of the reconstruction algorithm with respect t o a w hen the 
second line of (5.6) is deleted and (a, 6) — (a^°\0^) is fixed. Figure 6.3 'c) corre- 
sponds to the precise knowledge of the measurement setting, i.e., again ignoring the 
second line of (5.6), but fixing (a,0) = (/?,#) to be the parameter values describing 
the target configuration. Finally, the reconstruction in Figure 6.3 d) visualizes simul- 
taneous retrieval of the admittivity distribution and the measurement geometry by 
the whole two-stage algorithm of Section |5.3| 

The conclusions about the functionality of the different approaches are the same as 
in the previous experiment: The simultaneous retrieval of the admittivity distribution 
and the measurement geometry provides a reconstruction that is comparable to the 
case that the object shape and electrode locations are known accurately. On the 
other hand, ignoring the uncertainties in the measurement configuration gives a poor 
outcome. 

7. Concluding remarks. We have presented the Frechet derivative of the mea- 
surement map of practical EIT with respect to the (exterior) object boundary shape 
as a part of the solution to a certain elliptic boundary value problem. Through 
three-dimensional numerical studies based on simulated data, we have demonstrated 
that utilizing such a geometric derivative, the estimation of the object shape and 
the electrode locations can be incorporated into a Newton-type output least squares 
reconstruction algorithm in the framework on the CEM of EIT. 
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